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Abstract 

We consider formation of dissipationless shock waves in Bose-Einstein condensates 
with repulsive interaction between atoms. It is shown that big enough initial inho- 
mogeneity of density leads to wave breaking phenomenon followed by generation of a 
train of dark solitons. Analytical theory is confirmed by numerical simulations. 

Experiments on free expansion of Bose-Einstein condensate (BEC) have shown 1] that 
evolution of large and smooth distributions of BEC is described very well by hydrodynamic 
approximation 2] where dispersion and dissipation effects are neglected. At the same time, 
it is well known from classical compressible gas dynamics (see, e.g. jS]) that typical initial 
distributions of density and velocity can lead to wave breaking phenomenon when formal 
solution of hydrodynamical equations becomes multiple-valued. It means that near the 
wave breaking point one cannot neglect dispersion and/or dissipation effects which prevent 
formation of a multiple-valued region of a solution. If the dissipation effects dominate the 
dispersion ones, then the multiple- valued region is replaced by the classical shock, i.e. narrow 
layer with strong dissipation within, which separates smooth regions with different values of 
density, fluid velocity and other physical parameters. This situation was studied in classical 
gas dynamics and found many practical applications. If, however, the dispersion effects 
dominate dissipation ones, then the region of strong oscillations is generated in vicinity of 
wave breaking point jUE|. Observation of dark solitons in BEC 00 El shows that the main 
role in dynamics of BEC is played by dispersion and nonlinear effects taken into account by 
standard Gross-Pitaevskii (GP) equation |5J, and dissipation effects are relatively small and 
can be considered as perturbation. Hence, there are initial distributions of BEC which can 
lead to formation of dissipationless shock waves. Here we shall consider such possibility. 

The starting point of our consideration is the fact that the sound velocity in BEC is 
proportional to the square root from its density (see, e.g. [5] and references therein). Thus, 
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if we create inhomogeneous BEC with high density hump (with density ~ p\) in the center 
of lower density distribution (with density ~ po), and after that release this central part of 
BEC, then the high density hump will tend to expand with velocity ~ ^fpi greater than 
the sound velocity ~ v /po of propagation of disturbance in lower density BEC. As a result, 
wave breaking and formation of dissipationless shock wave can occur in this case. Note that 
initial distribution of this type was realized in the recent experiment JU] where generation 
of oscillations was also observed. In this experiment the lower density disk-shaped BEC 
was confined by magnetic trap and density distribution had standard Thomas-Fermi (TF) 
parabolic form. Additional potential was applied to BEC in the central part of above TF 
profile which lead to narrow parabolic hump in the density distribution. After the central 
potential was switched off, the hump started to expand against wide lower density TF 
profile leading to generation of oscillations in the transition region between high and low 
density condensates. The theory of dissipationless shock waves in media described by one- 
dimensional (ID) nonlinear Schrodinger (NLS) equation was developed in [HUE]- Since the 
GP equation in some cases can be reduced to the ID NLS equation, this theory can be applied 
to description of dissipationless shock waves in BEC. Here we suggest such description and 
confirm it by numerical simulations. 

We consider BEC confined in a "pancake" trap with the axial frequency uj z much greater 
than the transverse one u x — u y — Uj_. Let the density of atoms in the central part of 
BEC be of order of magnitude n and satisfy the condition riQa s a 2 -C 1, where a s > is the 
s-wave scattering length and a z = {h/ma z ) l l 2 is the amplitude of quantum oscillations in 
the axial trap. Then the condensate wavefunction ip can be factorized as ip — <p(z)ty(x,y), 
where <j)(z) = ix^^a^ 1 ^ 2 exp(— z 2 /2a 2 z ) is the ground state wavefunction of axial motion, and 
\l/(x, y,t) satisfies the reduced 2D GP equation 

h 2 

ifr% = ~y~ + Vw) + V ( x i V)* + 92d\^\ 2 ^, (1) 

where V(x,y) is the potential of a transverse trap, #2D = 2\/2Tih 2 a s / {ma z ) is the effective 
nonlinear interaction constant, and \l/ is normalized to the number of atoms, J \^>\ 2 dxdy = N. 
The initial distribution of density is determined by the potential V(x, y) and consists of 
wide background with a hump in its center. We assume that the background width is much 
greater than the hump's width, so that at the initial stage of evolution we can consider an 
expansion of the central part against the constant background. In a similar way, at the 
initial stage of evolution, when the radius of the central part does not change considerably, 
we can neglect the curvature of axially symmetrical distribution and consider its ID cross 
section. It means that we can neglect the dependence of \1/ on y coordinate and consider \l/ 
as a function of x and t only. As a result, we arrive at ID NLS equation with inhomogeneous 
initial distribution of density. To simplify the notation, we introduce dimensionless variables 
t' = tu z t/2, x' = x/a z , u = (2v^2™ s a 2/ Mo) 1 / 2 \I/ . Then the initial stage of evolution of the 
wavefunction profile in the x axis cross section is governed by the NLS equation, 

iut + u xx — 2\u\ 2 u = 0, (2) 

where the primes in t' and x' are omitted for convenience of the notation. 
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Evolution of smooth pulses before the wave breaking point can be described in the hy- 
drodynamic approximation which can be achieved by substitution 



u(x, t) = \J p(x, t) exp ( i / v(x',t)dx' ) (3) 



and separation of the real and imaginary parts. As a result we obtain the system 

\p t + {pv) x = 0, \v t + vv x + p x = 0, (4) 

where we have neglected the so-called "quantum pressure" term with higher space derivatives 
what is correct until the density distribution has smooth enough profile. To get simple 
qualitative picture of the wave breaking of BEC density distribution, let us consider idealized 
case with a box-like hump in the initial distribution, 

p(s,0) = ("» ! X jj a ' (5) 
Pi, \x\ < a. 

Although this distribution has a parameter with dimension of length — width of the hump — it 
does not play any role for t < a/(2y / pi), i.e. until two waves propagating inward the hump 
with sound velocity C\ = lyfpi meet at x — 0. Hence, for this initial period of evolution the 
solution for x > can only depend on the self-similar variable £ = (x — a)/t, centered at 
x = a, and for x < on £ = (x + a)/t centered at x — —a. Since the picture is symmetrical, 
it is enough to consider only a half of the solution corresponding to x > 0. It is easy to find 
that the density is given by the formulae 

for < x < a — 2^/pTt, 
p r xt \) i«vhi a)/2t] 2 j9 for a - 2^/plt < x < a + 2( v /p7 - y/pfit, ^ 



)/2t] /9 for a + 2^p- t < x < 2 + 2( v ^7- y/pfit, 
for x > a + 2y / p^t, 

and similar formulae can be written for the velocity field v(x,t). This solution describes 
the wave breaking phenomenon which takes place if pi > 4p . The density profile shown in 
Fig. 1 clearly illustrates the origin of the multiple-valued region which should be replaced 
by the oscillatory shock wave when the dispersion effects are taken into account. 

For more realistic initial pulses the density profile is a smooth function without cusp 
points. In vicinity of the wave breaking point the solution can be approximated by a cubic 
function for one Riemann invariant A + = v/2 + J~p of the system (J3J) and by constant value 
for another one A_ —v/2 — J~p (see [TT] I12j). After Galileo and scaling transformations the 
hydrodynamic solution can be written in the form 

x - (3A+ + A_)t = -A+, A_ = const, (7) 

and again for t > it has a multiple- valued region of A+. It means that dispersion effects 
have to be taken into account which lead to formation of dissipationless shock wave after 
wave breaking point. 

In framework of Whitham theory of modulations ^3 IS] one can obtain an approximate 
solution of the NLS equation (J2J in analytic form where the dissipationless shock wave is 
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presented as a modulated periodic nonlinear wave solution of the NLS equation. The density 
is expressed in terms of Jacobi elliptic function 



p(x,t) = \u{x,t)\ 2 = i(Ai - A 2 - A 3 + A 4 ) 2 + (Ai - A 2 )(A 3 - A 4 )sn 2 [ v /(A 1 - A 3 )(A 2 - A 4 )f , m], 

(8) 

where 

c (\ \ \ i \ \+ (Ai — A 2 )(A 3 — A 4 ) 

f = x - (Ai + A 2 + A 3 + A 4 )t, m = — — — — , (9) 

[M — — A 4 J 

and parameters Aj, i = 1,2, 3, 4, change slowly along the dissipationless shock. Their depen- 
dence on x and t is determined implicitly by the solution 

x — Vi(X)t = Wi(X), i = l,2,3; A 4 = A = const (10) 

of Whitham equations, where Whitham velocities V{ and Wi are given by quite complicated 
expressions in terms of elliptic integrals (see [TT] I12j): 

Wi = -l 5 w® + |A^ (2) - + ^A 3 , i = 1, 2, 3, (11) 

w f = WW + (v i -s 1 )d i WV e \ (12) 
W W = V = Su W M = | s 2 _ i S2) ^(3) = a s 3 _ | SiS2 + i S3) (13) 

v i( X ) = ( 1 - 7TF 9 i V. di^d/dXi, i = 1,2,3,4, (14) 



where 

KM 



, = (15) 

V / (Ai-A 3 )(A 2 -A 4 ) 

is a wavelength, K(m) is the complete elliptic integral of the first kind, and s\, s 2 , s 3 are 
determined by the expressions 

s% = ^ Xj, s 2 = ^ XjXj, s 3 = XiXjX k . (16) 

i «<j i<j<k 

Equations (|TT|) can be solved with respect to Aj, i = 1,2,3, giving them as functions of x 
and t. Subsequent substitution of these functions A«(x, t), i = 1,2,3, into (jEJ) yields the 
modulated periodic wave which represents the dissipationless shock wave. The resulting 
profile of density in dissipationless shock wave is shown in Fig. 2. At one its edge it consists 
of the train of dark solitons, and at another edge describes small amplitude oscillations 
propagating with local sound velocity into unperturbed region described by smooth solution 
of hydrodynamical equations. The modulated periodic wave replaces the multiple- valued 
region shown by dashed line which was obtained in hydrodynamic approximation after the 
wave breaking point. This multiple- valued region is analogous to that in Fig. 1 with account 
of change of variables due to Galileo and scaling transformations. 
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To check the described above picture of formation of dissipationless shock wave and to 
extend it to real 2D situation, we have solved numerically the 2D generalization of Eq. (J2J), 



iu t + u xx + u yy — 2\u\ 2 u = 0, 



(17) 



with the initial condition 



p(r, 0) 



Po + (pi-po)(l-r 2 /a 2 ), 
Pa, 



r\ < a, 
r\ > a, 



(18) 



(where r = (x 2 + y 2 ) 1 ^ 2 ) similar to one studied experimentally ^Uj. Plot of two-dimensional 
initial density distribution is shown in Fig. 3 and density distribution after time evolution t = 
2 is shown in Fig. 4. As we see, the parabolic hump expands with formation of dissipationless 
shock wave in the transition layer between high density region to low density one. To see 
more clearly the evolution of the hump, its cross section profiles are shown in Fig. 5 at 
different values of time t. Slowly propagating dark solitons are clearly seen as well as small 
amplitude sound waves propagating into undisturbed low density region. Dissipationless 
shock wave generated at the right side of the profile coincides qualitatively with results of 
analytic theory shown in Fig. 2. One may assume that oscillations in BEC density profile 
observed in experiment ^0] have the same origin. 

Further evolution of the hump will ultimately lead to its spreading over large area with 
small amplitude oscillations, that is the shock wave does not persist permanently. Rather, 
it is a transient phenomenon caused by different values of characteristic velocities in high 
density hump and low density background. Slowly propagating dark solitons are generated 
in the transient layer between these two regions in order to reconcile two different values of 
velocities of disturbance propagation. This mechanism of dark soliton generation is quite 
general and can manifest itself in various geometries and initial BEC distributions l . 

In conclusion, we have studied theoretically and numerically the process of formation of 
dissipationless shock waves in the density distribution of BEC. We believe that oscillations 
in BEC density profile observed in recent experiment ^0] can be explained by this theory. 
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Figures captions 



Fig. 1. Wave breaking of a box-like initial density distribution (shown by dashed line) in 
hydrodynamic approximation. The plot corresponds to the solution (jUJ) with a = 10, po — 1j 
pi — 10, t — 1. Point A propagates inward the box with local sound velocity 2y/p~[, point C 
propagates outward along the background density with local sound velocity 2 V //^, and point 
B corresponds to intersection of two simple wave solutions with profiles AB and CB. 

Fig. 2. Formation of dissipationless shock wave after wave breaking point according 
to Whitham modulation theory applied to ID NLS equation. Dashed line corresponds to 
multiple-valued region arising in hydrodynamic approximation given by Eq. ((7|) (it is analo- 
gous to the region between the points B and C in Fig. 1), and solid line represents modulated 
periodic wave given by Eqs. (jBJ) and fTUJ) . Both profiles are calculated for A = —10 at time 
t = 1. 

Fig. 3. Two-dimensional initial distribution of BEC density with paraboloid hump on 
constant background given by Eq. (fTHj) with a = 10, po = 1, Pi = 10. 

Fig. 4. Two-dimensional density distribution of BEC after time t = 2 of evolution from 
initial paraboloid density on constant background according to numerical solution of 2D NLS 
equation. 

Fig. 5. Cross sections of density profile at different evolution time according to numerical 
solution of 2D NLS equation (fT7j) with initial condition (f*TS|) . 
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